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Disks moving in a narrow channel have many features in common with the glassy behavior of 
hard spheres in three dimensions. In this paper we study the caging behavior of the disks which 
sets in at characteristic packing fraction (pd- Four-point overlap functions similar to those studied 
when investigating dynamical heterogeneities have been determined from event driven molecular 
dynamics simulations and the time dependent dynamical length scale has been extracted from 
them. The dynamical length scale increases with time and, on the equilibration time scale, it is 
proportional to the static length scale associated with the zigzag ordering in the system, which 
grows rapidly above (pd- The structural features responsible for the onset of caging and the glassy 
behavior are easy to identify as they show up in the structure factor, which we have determined 
exactly from the transfer matrix approach. 

PACS numbers: 64.70.Q-, 05.20.-y, 61.43.Fs 


I. INTRODUCTION 

It has been frequently observed that disks moving in 
a narrow channel can provide useful insights into glassy 
behavior [1-7]. In a recent paper [6] two of us studied 
the static and dynamic properties of N disks of diam¬ 
eter CT, which move in a narrow channel consisting of 
two impenetrable walls (lines) spaced by a distance Hd, 
such that 1 < Hd/u <1-1- •\/5/2 (see Fig. 1). With 



FIG. 1: (Color online) Illustration of the disk-channel system 
as in Ref. [6]. The shaded disks represent a defect in the 
developing zigzag order. 


channels of this width only nearest-neighbor disks can 
interact with one another and neighboring disks can¬ 
not pass each other so that their ordering is preserved: 
0 < xi < X 2 < ■ ■ ■ < Xn < L, where Xi is the position of 
the center of disk i, measured along the channel, and L is 
the total length available to the disk centers. The coordi¬ 
nate Ui of the zth disk is measured from the center of the 
channel and can vary between ±/i/2, where h = Hd — cr. 
The packing fraction cp in our system of disks is defined 
as 


Nna'^ 


4HdL ■ 


( 1 ) 


Our disks move ballistically, colliding elastically with 
each other and the channel walls. It is useful for many 
purposes to regard the disks as being contained within a 
system of average length L by pistons at the ends of the 
channel which exert a force / to counteract the momen¬ 
tum transferred by the disks colliding with them. For 


any finite value of this force, large fluctuations in the x- 
coordinates of the disks cause the time-averaged density 
of disks to be independent of x, so that the system is 
never crystalline. The static, equilibrium properties of 
this simple system can be determined exactly by use of 
the transfer matrix [6, 8-11], but the chief purpose of this 
paper to discuss the dynamics. The dynamical properties 
of our system must be determined from simulations, and 
to this end we have used event-driven molecular dynamics 
to handle the collisions of the disks with each other and 
with the channel walls. We find interesting similarities 
and also some instructive differences with the dynamics 
of three dimensional glass-forming liquids. 

It was found some time ago that in the system of disks 
in a channel there is a packing fraction cpd ~ 0.48 above 
which the dynamics are activated [1]. Similarly, for hard 
spheres in three dimensions, the relaxation time grows 
rapidly above a packing fraction of « 0.58. For a re¬ 
view of this and other features of hard sphere systems see 
Ref. [12]. We have also observed that for (p > (pd ^ 0.48 
zigzag ordering of the disks starts to grow rapidly [6] ; that 
is, the onset of the slow dynamics is connected with the 
growth of this particular kind of order. There is a long 
tradition of associating glassy behavior with geometri¬ 
cal features associated with the arrangements of parti¬ 
cles around a given particle, see e.g. [13-17], and our 
work on disks in a narrow channel is entirely consistent 
with these ideas. In our earlier work [6, 7] a length scale 
^ associated with zigzag order has been determined from 
the decay with s of the correlation function 

iVi Vi+s) (-1)'' exp(-s/^). (2) 

It is a measure of the number of disks over which the 
zigzag order, a form of bond-orientational order, per¬ 
sists; it is of the same order as the separation of defects 
in the zigzag order, where the defects are of the kind 
illustrated in the top and bottom panels of Fig. 2. In 
three dimensional systems the ordering associated with 
glassy behavior is complicated [17]. Furthermore, the or¬ 
dering associated with glassy behavior is not apparent in 
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FIG. 2: (Color online) The transition state for motion of a 
defect. In the top and bottom diagrams, the two blue-shaded 
disks are a defect in the zigzag arrangement of disks. The de¬ 
fect can move when one disk crosses the channel by squeezing 
between its neighbors: the system passes through the tran¬ 
sition state shown in the middle diagram to reach the defect 
state shown in the bottom diagram. In the top diagram the 
defect involves disks 3 and 4; in the bottom diagram the de¬ 
fect involves disks 4 and 5, when the disks are numbered from 
the left. The net motion of the defect is to the right, and A;, 
is the extra length needed to allow this motion. 

changes to the structure factor as it is cooled through 
the glass transition. In Sec. II we calculate the structure 
factor of our system, defined as 

S{k^,ky) = ^'^{exp[ik^{xi - Xj) + iky{y, - yj)]) , (3) 
hi 

at packing fractions close to (fid by use of the transfer ma¬ 
trix formalism and so it is exact to numerical accuracy. 
It is found to change rapidly near <f>d- This marked dif¬ 
ference with the behavior in three dimensions where the 
structure factor hardly alters near the glass transition 
indicates that glass behavior in three dimensions must, 
as suggested in Refs. [13-17], involve higher order corre¬ 
lations which have little impact on the structure factor, 
unlike the simple zigzag orientational order of the narrow 
channel system. 

In Sec. Ill we shall study the onset of the slow dy¬ 
namics which sets in around (f>d- In three dimensions 


the slow dynamics are normally attributed to the on¬ 
set of caging behavior where a particle is trapped by its 
neighbors and this behavior is believed to be captured 
by the mode-coupling approach. The particle can escape 
its cage on the alpha relaxation time r^. To determine 
whether or not caging occurs in our system we have stud¬ 
ied the mean-square displacements 

1 ^ 

' i=l 

where the average is over different initial states. We find 
that there is caging of particles in our system, i.e. there 
is (at large enough packing fractions) a plateau in A^(t) 
before its long time limit is reached. We hnd that in our 
system there are two distinct large time scales, which we 
call T and td following Ref. [6]; we can obtain both of 
them from the behavior of A^(f). The smaller of these 
time scales, r, is the typical time scale for a disk to cross 
from one side of the system to the other by the process 
shown in Fig. 2. It marks the time at which the particle 
starts to escape from its cage or the end of the plateau in 
AHt): in three dimensional systems this would be called 
the alpha relaxation time . At packing fractions above 
(f>d, escape from a cage can be regarded as a displace¬ 
ment of a defect in the zigzag ordering of the disks as is 
also illustrated in Fig. 2. We are able to determine from 
the simulation the diffusion constant for the movement 
of such defects. Its dependence on packing fraction is 
consistent with the general picture of a crossover from 
fragile glass behavior at low packing fractions to strong 
glass behavior at high packing fraction, which has been 
investigated previously in Refs. [2, 4]. 

The second long time scale tjj is essentially the longest 
time scale in the system. In equilibrium defects are ther¬ 
mally nucleated in pairs and the defects produced then 
diffuse and annihilate with each other. It is this process 
of diffusion with creation and annihilation which takes 
place on the time scale T]j . There is a simple relation be¬ 
tween T and td: td ~ [6]. Note that ^ is the typical 

distance between defects in equilibrium, so that /D is 
the time it would take for a defect to move a distance 
this gives a diffusion coefficient D for the defects which 
varies as 1 /r [6] . T]j is determined from the time at which 
A^(f) approaches its equilibrium limit. 

A characteristic of glassy dynamics is the appearance 
of a plateau in the decay of certain time-dependent corre¬ 
lation functions. This plateau eventually decays to zero 
after the time the alpha relaxation time. The time to 
reach the plateau is the beta relaxation time Correla¬ 
tion functions with this feature have not previously been 
studied in our narrow channel system. The existence of a 
plateau in A^(f) implies that there will be such a plateau 
in the decay of 

Rit) = ^ Viit )); (5) 

i 

at long times R{t) approaches zero. As noted above. 
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our results support a model in which the diffusion of de¬ 
fects leads to equilibrium. In models of this kind, R{t) 
is expected to decay with a stretched exponential form 
exp(—where td is the equilibration time of 
the system [18]. 

We have also studied a correlation function related to 
the X 4 (t) which has been much studied in three dimen¬ 
sions [15, 19, 20]. In three dimensions rises to a 
peak at the alpha relaxation time and then decays to 
zero as the particles escape their cages. Our X 4 (t) in¬ 
creases to a constant, plateau value as t increases and 
stays at this value as t —?► oo. In this regard the nar¬ 
row channel system behaves more like a spin glass than 
a structural glass [21]. We have also extracted a dynami¬ 
cal correlation length ^^{t) from a four-point correlation 
function. In our system as t —)■ td , ^4 grows towards the 
static correlation length which measures the extent of 
zigzag, i.e. structural, order. 

In Sec. II we show how the structure factor can be 
determined exactly from the transfer matrix procedure 
and we demonstrate that it changes rapidly for densities 
close to (j)ci. Our dynamical studies are in Sec. III. We 
conclude with a discussion in Sec. IV. 
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FIG. 3: (Color online) The structure factor S{kx,ky = 0) 
as a function of Pfa and for the narrow-channel system 
with h = y/Zal2. The position of the first peak in S (near 
kx<y = 6 for Pfa = 4) changes rapidly for Pfa in the range 
6.5 to 8, following the change in periodicity that accompanies 
the growth of zigzag order. Adjacent contours are spaced by 
0.1 in logjQ S and contour lines with S' = 1 are marked in red. 


II. STRUCTURE FACTOR 

The structure factor is well known from the part it 
plays in the scattering of electromagnetic radiation by 
liquids [22]. It is the Fourier transform of the density- 
density correlation function and so provides information 



FIG. 4: (Color online) The structure factor S{kx,ky = ir/h) 
as a function of pPa and k^cr for the narrow-channel system 
with h = y/3af2. A peak near k^a = 4 grows rapidly for 
Pfa > 4, consistent with the growth of zigzag correlations. 
Contour lines are described in the caption to Fig. 3. 

on the relative positions of scatterers within the system. 
Like a liquid, our system of disks in a channel has no 
long-range order (along the channel) for finite values of 
the longitudinal force /; but, unlike in the bulk of a liq¬ 
uid, the short-range order is strongly affected by the pres¬ 
ence of confining walls, leading to the zigzag correlations 
discussed in Sec. I. In this section we show that the struc¬ 
ture factor can be calculated essentially exactly for disks 
in a narrow channel. Our numerical results for the case 
h = •\/3o’/ 2 show a rapid change in the short-range or¬ 
der for Pfa in the range 6.5 to 8 , which corresponds to 
packing fractions p in the range 0.45 to 0.50. This range 
of packing fractions correlates closely with the onset of 
activated dynamics, as we discuss later in Sec. III. 

In the limit V —k 00 , the definition of the structure 
factor, Eq. (3), may be rewritten in the form 

00 00 

S{ky,,ky)^ Sr^ = l + 2ReJ2Sn, ( 6 ) 

n— — (x> n—1 

where 

_ ^^iik^[x„-xo]+ky[y„-yo])\^ ^ 

At ky. = 0, S{kx, ky) has a delta-function singularity and 
the sum on the right-hand side of ( 6 ) diverges. As we now 
show, for ky ^ 0 the sum can be evaluated relatively sim¬ 
ply by solving a pair of integral equations. One of these 
equations is known from the transfer-matrix formalism 
introduced by Barker [ 8 ] and applied by Kofke and Post 
[9] to the problem of hard disks in a channel. We follow 
the latter authors in using an ensemble in which the lon¬ 
gitudinal force / is constant and we refer the reader to 
their paper [9] for details. 
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Let tpn{y) be the eigenfunctions of Kofke and Post’s 
integral equation 

fh/2 

A„V’n(yi)= / (8) 

J-h/2 

where yo and yi are the y-coordinates of a neighboring 
pair of disks and ai^ — [a^ — {yi — yo^Y^^ is the dis¬ 
tance of closest approach of their centers, measured along 
the cc-axis. Approximations to the eigenfunctions ipn and 
the eigenvalues A„ can be found by discretizing Eq. (8), 
which converts it to a real-symmetric matrix eigenvalue 
problem. The eigenfunctions (taken to be real) can be 
normalized so that 

j[>Pn{yi)f dyi = I. (9) 

In this and subsequent equations, the limits of the y in¬ 


tegration are —h/2 and h/2, the same as in Eq. (8). 

Equilibrium expectation values, such as those needed 
for the quantities Sn defined in (7), can be expressed 
as integrals involving the eigenfunction tpi which corre¬ 
sponds to the largest eigenvalue Ai. We illustrate this for 
the calculation of below. 

5'i is the expectation value of exp{ikx [a;i —a^o] +iky [yi — 
?/o]), which is a function of yo, yi and the gap s between 
disks 0 and 1, defined by 

S + (Tifi = Xi - Xq . (10) 

From the results of Ref. [9] , the probability distribution 
for the variables yoi 2 / 1 ; and s is proportional to 

■*/'i(2/i)e"'^-^^®+‘^"-‘’' V’i(2/o) ■ 


Accordingly, S'! is given by 


Si = 


^ik^(s+(7i^o)+iky(yi-yo 




I Myi) I /o°° Myo) ds dyp dyi 

/ fpiiyi) f ipiiyo) ds dyo dyi 


( 11 ) 


After completing the integrals with respect to s and using the eigenvalue equation (8) and the normalization condition 
(9) to simplify the denominator, we obtain 


Si=J Myi) { x,{pf-tk,) J 

= J tpiSipidyi, (12) 


in which the bracketed expression in the first line defines 
the action of the integral operator S on -tpi. More gener¬ 
ally, for n > 1 one can write 

Sn = dyi , (13) 

so that the sum in Eq. (6) becomes 

OO OO p n 

'^Sn = '^ / i’l S'^'fpidyi = / ipiScjidyi, (14) 

where cj) is the solution of 

(j) = lj;i+ S(j), (15) 

which is a Fredholm equation of the second kind. Given 
the function ipi found by solving the discretized Eq. (8), 
the calculation of (j) requires only the solution of the set 
of linear equations obtained by discretizing Eq. (15). 

Finally, in terms of 'ipi{yi) and (piyi)-, the structure 
factor is given by 

J ipi S(j)dyi , (16) 


which depends on kx and ky via (j) and the operator S. 

As indicated above, the equations (8) and (15) can be 
solved by discretization. For all but the smallest values 
of Pfa, the function tpi is concentrated near the walls 
and it is sampled at smaller intervals in those regions, to 
keep the dimension of the matrices relatively small. To 
implement this nonuniform sampling we make the change 
of variable 


^ sinh(a/3/y) 
sinh(a/3/h/2) ’ 


(17) 


where a = \hl\J and the values of u are taken to 
be uniformly spaced in the interval [—1,1]. The symme¬ 
try of the transfer matrix is preserved by solving Eq. (8) 
for the function rather than i/ji. Matrices 

of dimension 100 x 100 are sufficient to reproduce the 
results presented in this Section. 

Numerical results are shown in Figs. 3-7. In Fig. 3, 
the structure factor is plotted as a function of /3/cr and 
kx for the case ky = 0. For ky = 0, the structure factor 
is sensitive only to correlations in the j/-averaged density. 
Zigzag order is growing rapidly for f3fa in the range 6.5 


S{kx,ky) = l + 2Re 
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to 8, and this can be seen in Fig. 3 as a rapid change in 
the nature of the first maximum with respect to kx- For 
very small values of /, the oscillations in S{kx,ky = 0) 
are small and (much as for a low-density gas of hard 
rods) their presence is due to the discontinuity of the pair 
distribution function g{r) at r = cr. But as / increases, 
a new peak emerges and rapidly becomes dominant. Its 
position reflects the changing spatial periodicity in the x- 
direction: for /3fa —>■ oo, the peak approaches kx = 47r/cr, 
which can be understood as 27r/(tT/2), where (t/ 2 is the 
spatial period of the y-averaged density in this limit. 

Figure 4 shows the evolution of the structure factor for 
the case ky = ir/h, where the value of ky has been cho¬ 
sen to better illustrate the growth of zigzag order. Note 
that when Pfa is large, nearest-neighbor disks are sepa¬ 
rated by Aa; a 12 along the channel and by Ay « ±/i 
in the transverse direction. For ky = tt/Ii, the nearest 
neighbors will scatter in phase when kx Ax « tt, i.e. for 
kx ~ 21:la. This scattering results in a peak in the struc¬ 
ture factor that is clearly visible in Fig. 4. The peak first 
appears near kxa = 4 for Pfa « 4. It strengthens as Pfa 
increases and it eventually approaches kxa = 27r. 
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FIG. 5: (Color online) The structnre factor S{kx,ky) for 
Pfa = 5 for the narrow-channel system with h = \fZal2. 
Contour lines are described in the caption to Fig. 3. 


Figures 5-7 show the structure factor for three val¬ 
ues of the force for which S{kx,ky) is changing most 
rapidly: Pfa = 5, 7.5, and 10. The maxima in S grow 
and become narrower (in the kx direction) as / increases 
and the zigzag correlations strengthen. For very large /, 
the widths of the peaks decrease as {Pfa)~’^. We note 
that this /-dependence of the peak-width is the same as 
is found for a one-dimensional gas of hard rods, whose 
structure factor was derived analytically by Zernike and 
Prins [23]. 



FIG. 6: (Color online) The structnre factor S{kx,ky) for 
Pfa = 7.5 for the narrow-channel system with h = y/ial2. 
Contour lines are described in the caption to Fig. 3. 
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FIG. 7: (Color online) The structnre factor S{kx,ky) for 
Pfa = 10 for the narrow-channel system with h = y/Zal2. 
Contour lines are described in the caption to Fig. 3. 

III. TIME-DEPENDENT BEHAVIOR 

In order to study time dependent effects we have used 
event driven molecular dynamics based upon the code 
referred to in Ref. [24]; the speed of the program was 
improved by using the fact that in our narrow-channel 
system only nearest-neighbor disks can collide. The ini¬ 
tial state from which the system evolves was created 
by means of the Lubachevsky-Stillinger algorithm [25], 
starting from a random configuration of small disks. 
Their diameters were slowly increased to the desired 
value during the course of a simulation which preceded 
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the long runs used to study dynamics in the equilibrated 
system. 

The force / along the channel was computed by using 
a virial formula. Suppose that at the instant of collision 
between two disks, the ^-separation of their centers is 
Ax(c) > 0, where c labels the collision. The x-component 
of the momentum transferred from the disk at smaller x 
to the disk at larger cc is a positive quantity, Apx{c) > 
0. With the values of Ax{c) and Apx{c) determined by 
simulation, the average longitudinal force can be found 
from 

fL=^ +—^Ax{c)Ap^{c), (18) 

P "^sim ^ 

where the sum includes all disk-disk collisions c that oc¬ 
cur during the simulation time Tgim- 

All of our simulations were performed with N = 10000, 
h = and /3 = cr = m = 1, where m is the mass of a 

disk. To check their accuracy we calculated the equation 
of state and compared it with the results of the exact 
transfer matrix calculation [6], as shown in Fig. 8. The 



FIG. 8: (Color online) Equation of state Pfa versus ij) for 
the narrow-channel system with h — \/3(t/ 2. The solid line 
shows the equation of state calculated from the transfer ma¬ 
trix results of Ref. [6], while the data points were obtained 
from long-time simulations of A = 10000 disks. 

agreement is excellent up to <() = 0.65. There is a dis¬ 
crepancy at (p = 0.70 (but too small to be visible in the 
figure) as at this packing fraction the time scale for equi¬ 
libration, T£), is becoming comparable to our simulation 
time. 


A. Time scales 


To see the emergence of caging behavior it is conve¬ 
nient to study 


A^{t) 



2=1 


1 


(19) 


where the angled brackets indicate an average over runs 
with different initial states. Though it appears unnatural 
at first sight, this quantity was introduced in Ref. [26] 
to miminize the contribution of rattlers. The ability of 
rattlers to move a distance of O(tT) tends to dominate 
the mean-squared displacements at small times t. We 
have a similar problem here in that a few disks which 
border gaps are able to cross from one side of the channel 
without hindrance. In Fig. 9 one can see the emergence of 
a plateau as p increases from 0.45 to 0.50, which suggests 
that caging of the disks is setting in as the zigzag order 
develops around (pd- It was in Ref. [1] that (pd ~ 0.48 was 
first identified as the onset point for activated dynamics. 

It is easier to understand the behavior of the un¬ 
modified mean-square displacement A^(t) as defined in 
Eq. (4). However, Fig. 10 shows that for this quantity the 
plateau is only clearly visible at values of <p above 0.60, 
which is well into the density regime where the dynamics 
are activated [1, 6]. 



FIG. 9: (Golor online) Time dependence of the mean-squared 
displacements (in units of cr^), of 10000 disks defined via 
Eq. (19), a quantity which was introduced in Ref. [26] to mini¬ 
mize the contribution from the most mobile disks. The results 
for the various packing fractions p which were studied are av¬ 
eraged over 20 independent equilibrated trajectories. Time is 
in units of . 

Some of the features on display in Fig. 10 include: 

1. Final long-time limit: In the limit t —)■ oo, A^(t) 
reaches a finite value. (In three dimensions in the 
same limit, the mean-square displacement of a par¬ 
ticle increases without limit linearly with t.) From 
its definition, A^ = {yf{t)+y‘p{0) — 2yi{t)yi{0))-, the 
last term in the brackets gives 2R{t), which tends 
to zero in the long time limit, so A^ tends to 2{yf), 
a quantity which is close to 2(/i/2)^ in the limit of 
large / when the disks are mostly pushed to the 
sides of the channel. 

2. For small t, A^(t) increases as i.e., the motion 
is ballistic. This regime is larger at smaller values 
of p, for which the gaps between disks are larger. 











7 



logic t 

FIG. 10: (Color online) Time dependence of the mean- 
squared displacements (in units of cr^), for trajectories of 
10000 disks at varying packing fractions (j). Time is in units 
of The dashed lines fit the central regions where 

the gradient is approximately 0.5, indicating a process of re¬ 
laxation that is dominated by the diffusion of defects, as 
suggested in [6] as a mechanism for the a-relaxation. For 
<!> ^ 0.60 a plateau is seen to form between the small-t ballis¬ 
tic and large-t diffusive regimes. The plateau corresponds to 
disks becoming trapped in cages formed by their neighboring 
disks, which require cooperative motion to break. At very 
long times A^(t) tends to a constant value, 2{yf). 

3. A “shoulder” begins to form around (p ~ 0.60, and 
develops into a clearly visible plateau for (p > 0.65. 
This is a clear analogue of the caging effect seen in 
three dimensions. It sets in at higher packing frac¬ 
tions than (pd ~ 0.48, the packing fraction at which 
the growing zigzag order results in the dynamics 
becoming activated. 

4. Beyond the “glassy” plateau for (p > 0.60, and 
above the ballistic regime for cp < 0.60, there is 
a time scale r beyond which a power-law depen¬ 
dence on t sets in, A^(t) oc This dependence 
on t is due to the slow diffusion of defects in the 
zigzag arrangement of disks. 

We can explain some of these features in greater detail. 
While activated dynamics may set in around a packing 
fraction p « 0.48, some disks still find at this density 
that they can easily cross the channel, and it is not until 
a packing fraction of 0.60 that the numbers of these “rat¬ 
tling” disks become negligible. The plateau represents a 
clear caging effect, and it lasts for a time t, the time it 
typically takes for a disk to cross from one side of the 
channel to the other by the transition state mechanism 
depicted in Fig. 2. Once the zigzag order sets in, the mo¬ 
tion of the defects as in Fig. 2 dominates the behavior of 
most of the disks. However, at packing fractions around 
0.48 there are still some disks that can travel from one 
side of the channel to the other with little hindrance from 
their neighbors. (A similar observation has been made in 
Ref. [4].) Their contribution to A^(t), defined in Eq. (19) 


is small, which enables one to see the emergence of the 
caging behavior in it at lower packing fractions than for 

The time scale r for a defect to move as shown in Fig. 2 
was studied numerically in Ref. [1] and explained using 
transition state theory in Refs. [6, 27]. At high packing 
fractions 

r ~ Toexp(/3/A{,), (20) 

where tq is of the order of a disk collision time. The ar¬ 
gument of the exponential in Eq. (20) can be understood 
from Fig. 2: PfAb is the work done against the piston in 
creating the extra length A;, in the system which allows 
the defect to move. In Ref. [6] it was shown that this ex¬ 
tra length was A;, = V4(7^ — — a —y/a — which can 
also be understood by inspection of Fig. 2. The plateau 
visible at larger packing fractions in Fig. 10 will end on 
the time scale r. 

We now turn to the details of the diffusive behavior, 
indicated by the dashed lines in Fig. 10. A^(t) increases 
as (yi(t)yi(O)) goes to zero. The crossing of disks from 
one side to the other of the channel is what drives this 
correlation function to zero. Fig. 2 shows that this hap¬ 
pens where there are defects in the zigzag arrangement 
of the disks. Let 6 be the concentration of defects, so 
that the number of defects is N6. This number is readily 
determined in numerical work using the procedure given 
in Ref. [3]. Figure 11 shows 0 as a function of packing 
fraction p obtained from the simulations and compared 
with the analytical approach of Ref. [6], which becomes 
exact at large packing fractions. As time increases the 
number of disks flipped by the diffusion of the defects 
will be of order N9p'~Di, where D is the diffusion coeffi- 
cent of a defect. It was argued in Ref. [6] that at large 
/, Dtd ~ 1/6^. Then 

t£, - To exp(/3/Ac), (21) 

where Ac = V+ cr — 3V The physical 
significance of the time scale td is that it is the time 
scale on which diffusing defects meet and annihilate [6]. 
It appears to be the longest time scale in the system - the 
time scale for full equilibration. Note that in the diffusive 
region 

A^ {t) ~ hP^ yJt/TD ■ (22) 

By extending the dashed lines in Fig. 10 to the points 
where they meet the analytic final values of A^(t) we 
can estimate values for td. Values of toItq are plot¬ 
ted in Fig. 12. The collision timescale tq was estimated 
from the mean collision rate per disk determined from 
our simulations, i.e. 

- = ^, (23) 

Tq I\ Tgini 

where is the total number of collisions that occurred 
during the simulation time Tsim- The agreement with 











FIG. 11: (Color online) Variation of the defect density 6 with 
packing fraction (j>. The points give the observed defect densi¬ 
ties in our simulations and the dashed line is the approximate 
theoretical relation predicted in Ref. [6] . The approximate re¬ 
lation is expected to improve in the limit of <(> — ^ <(>max, which 
is seen here in the increased agreement between simulation 
and analytic values as (j) increases, with the exception of the 
result for (j> = 0.7 where we found 4 times more defects than 
expected. This discrepancy is most likely due to poor equi¬ 
libration as the trajectory was simulated for a time period 
much shorter than its relaxation time. Relaxation times to 
are shown in Fig. 12. 


the prediction (21) is satisfactory for the results shown 
in Fig. 12. 

We have also used the results shown in Fig. 10 to ob¬ 
tain the diffusion coefficient for defects. At high den¬ 
sity, the diffusion coefficient should be related to r by 
Dtq ~ To/r, and so might be expected to provide confir¬ 
mation of Eq. (20). To find D, we make use of (22) in 
the form 

A^{t) ^ By/Di. (24) 

Extrapolation of the dashed lines in Fig. 10 back to t = 1 
gives, via (24), an estimate of 0y/~D. This in turn pro¬ 
vides D when we make use of the values of 9 found from 
our simulations. Results for Dtq obtained in this way 
are plotted in Fig. 13. The results show a significant de¬ 
parture from Eq. (20), but there are several reasons why 
we might expect our procedure to give poor results for 
the densities of interest here. First we note that for the 
smaller values of the packing fraction the lines of slope 
0.5 in Fig. 10 are not very convincing fits to the data: 
the linear portions of the curves are very short. Sec¬ 
ondly, the mean spacing of defects given by \/9 is not 
large (see Fig. 11) for (j) in the range 0.4 to 0.6: inter¬ 
actions between defects may well modify the diffusion 
coefficient significantly, leading to a 0-dependent factor 
in the relation D ^ 1/r. Finally, as shown in Ref. [6], 9 
itself is not a simple exponential function of (3 fa at these 
moderate values of (j), but instead can be calculated quite 
accurately (see Fig. 11) from a law of mass action. Thus, 


even if r has the activated form (20), it is not certain 
that this can be ascertained from our calculation of D at 
the moderate densities accessible via our simulation. 

Fortunately, it is not necessary to rely on an estimate of 
the diffusion coefficient to verify the activated behavior 
of T. Bowles and Saika-Voivod [1] have made a direct 
determination of the channel-crossing time from their 
molecular dynamics simulation. As shown in Ref. [6], 
their results are in satisfactory agreement with Eq. (20). 



FIG. 12: (Color online) Rapid increase in the ro-relaxation 
times with longitudinal force /, obtained by assigning a linear 
fit to the diffusive region of the mean square displacements 
(see Fig. 10) and extrapolating to find the time where it meets 
the analytic final values of A^(t). The dashed line fits the ex¬ 
ponential trend predicted by Eq. (21) that is In to (to ~ PfAc 
and becomes an increasingly better fit as / becomes large, i.e. 
in the limit (j) —>■ ())max. The circled data point, corresponding 
to (j) = 0.7, deviates from this trend: as noted in the caption 
to Fig. 11, the system failed to reach equilibrium for this value 
of (j). 

Finally we comment very briefly on the time-dependent 
correlation function R{t) = {yi(t) yi{0)), mentioned in 
Sec. I. R{t) decays to zero because of the diffusion of 
the defects in the regular zigzag arrangement of disks in 
the channel. Such a diffusive mechanism has been much 
studied [18] and leads to a stretched exponential decay 
R{t) ~ exp(—i/t/ru). We have not attempted a direct 
verification of this behavior of R{t), as it is expected only 
for very large times, t ^ to- 

B. Overlap correlations 

We turn now to a study of overlap correlations, which 
correlate the configuration of the system at time t with 
a configuration drawn from the equilibrium ensemble at 
time t = 0. Studies of such correlations reveal the exis¬ 
tence of dynamical heterogeneities in three-dimensional 
glass-forming liquids. It is therefore of interest to see 
whether the system of disks in a channel shows similar 
behavior. In particular we shall make use of the (self) 
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FIG. 13: (Color online) Variation of diffusion constant D oc 
1/r with longitudinal force / using the j/-intercept of linear 
fits to the diffusive regions shown in Fig. 10. The dashed line 
fits the exponential trend predicted by Eq. (20), Inr/ro ~ 
PfA,. 

overlap function 

(25) 

i 

where w[y^{t),y^{Q)\ = i(sign(y,(t)) sign(?/i(0)) + 1) is 
unity if disk i is on the same side of the channel at 
times 0 and t, and is zero otherwise; in the terminology 
of Ref. [20], this quantity is the mobility of disk i. Similar 
overlap functions have been studied in three dimensions: 
in that case, the overlap w[ri,rj] has been taken to be 1 
if |ri — Tjl < 0.3CT and zero otherwise [14, 15, 19]. A sig¬ 
nificant difference between this latter definition and our 
own is that our overlap function does not depend on the 
cc-coordinates of disks. This modification eliminates an 
effect, specific to our one-dimensional problem, of large 
fluctuations (^iV^/^) in the ^-coordinates of disks. These 
fluctuations can cause the overlap between two configu¬ 
rations to be small even in cases where the configurations 
would be identical, when described in terms of defects in 
the zigzag arrangement of disks. 

A quantity much studied for three-dimensional systems 
is the four-point susceptibility X4(^)) which is defined in 
terms of the variance of Q(t) via 

Xi{t)/N = {Q{tf) - {Q{t)f. (26) 

This has been calculated for our one-dimensional system 
and is shown in Fig. 14. In three dimensions it reaches a 
maximum on the timescale Tq, and subsequently decreases 
towards zero, but for our system there is no decay back 
to zero. This difference between dimensions c? = 3 and 
d = 1 has been discussed previously in Ref. [28] and can 
be understood qualitatively as follows. In three dimen¬ 
sions the self-overlap r(;[ri(t), ri(0)] becomes small and is 
likely to remain small once a sphere has escaped from 


its cage; it follows that the fluctuations in Q{t) will also 
be small for times that are long enough for most of the 
spheres to have escaped from their cages. This argument 
does not apply to our system of disks in a channel because 
a disk can cross the channel any number of times, so that 
w[yi(t), ?/i(0)] is a fluctuating quantity of order unity. For 
any value of t, the ^-coordinates of disks are correlated 
over a range ^ and are approximately independent over 
larger distances. From this we expect the contributions 
to Q{t) from independent regions of size ^ to be fluctu¬ 
ating quantities of order £^/N for times t > tjj. If the 
number of these regions is ~ we expect 

Xa/N = Var Q{t) ~ {^Nf , (27) 

which gives xa ^ P for long times. This result can be 
refined by using Eq. (2) to evaluate the right-hand side 
of (26); we obtain 

XA{t) « e/4 (28) 

for t ^ Tjj and ^ 1- The last estimate (28) is con¬ 

sistent with the results shown in Fig. 14 for (p = 0.49 
and 0.59, for which ^ = 3.9 and 20.9, respectively. 



Iogio t 


FIG. 14: (Color online) Xi{t), as defined by Eq. (26), as a 
function of time t for various packing fractions (f>. It ap¬ 
proaches its largest values at a time corresponding to td ■ 

We have also determined the dynamical length scale 
from the four point correlation function S 4 {kx, t) defined 
as 

SA{K^,t) = ^{Q{Kx,t) Q*{Kx,t)) , (29) 

where 

Q{i^x,t) =^e~^^^^yj{t)yj{Q) (30) 

d 

and Kx = 2'Km/N^ where m = 1, 2, ..., TV — 1. The 
Kx-dependence of SA{Kx,t) follows a roughly Lorentzian 
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form, as found in Refs. [15, 19], i.e., 




m 

l + nlilit) 


for Kx 7 ^ 0; the dynamical length scale ^ 4 (t) is obtained 
from a fit to the simulation data. 

We present results for ^ 4 (t) at two physically signif¬ 
icant values of f. At t = td the system has reached 
equilibrium, so that 

(31) 

p,<i 


Then, by using Eq. (2), one finds the correlation length 



0 2 4 6 8 10 




Utd) « e/2. (32) 

Here ^ is the number of disks for which zigzag order per¬ 
sists. This dimensionless quantity is calculated from the 
two largest eigenvalues, Ai and A 2 , of the integral equa¬ 
tion (8), using the expression 

e=l/ln(Ai/|A2|), (33) 

which follows from the spectral representation of correla¬ 
tion functions, discussed for the case of the Ising model 
in Sec. 2.2 of Ref. [29] and applied to a system of hard 
disks in a channel by Varga et al [10]. The results shown 
in Fig. 15 are in good agreement with the prediction (32). 

If we choose for t a value for less than the equilibration 
time td, then ^ 4 {t) will be less than ^/2. On the time 
scale T, the time for which particles are caged, 54(r) is 
not proportional to ^ and the results in Fig. 16 show it 
instead to be of order 1 and approximately independent 
of (f)toT (j) > (pd- This behavior of ^ 4 (r) is understandable, 
as on the time scale r the active regions are centered on 
the defects, which involve 0(1) disks. 

The above studies show that our system has dynam¬ 
ical heterogeneities of the kind expected in the defect- 
mediated scenario for glassy dynamics [28, 30]. In the 
packing-fraction regime where the dynamics are acti¬ 
vated, most of the disks will be largely frozen except 
those in the vicinity of a defect. The spacing between 
defects is of order ^ which can become very large. The 
defects move about and will eventually annihilate with 
each other. This happens on the time scale td, which 
is also the time scale on which new pairs of defects are 
typically nucleated [6]. 

IV. CONCLUSIONS 

There are some notable similarities between the behav¬ 
ior of our system of disks in a narrow channel and that 
found for a fluid of hard spheres in three dimensions. 
In each of these systems at higher densities the mean- 
squared displacement of particles has a plateau that per¬ 
sists up to the characteristic time for the breaking of 


FIG. 15: (Color online) Plot of the dynamical length scale 
ii(TD) versus the static length scale ^ of zigzag order. The 
gradient of the dashed line is 0.5, corresponding to the pre¬ 
dicted behavior, Eq. (32). 
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FIG. 16: (Color online) Plot of the dynamical length scale 
^ 4 (t) versus the packing fraction (f>. 


cages, r or Tq. Also, in both systems there is a crossover 
from non-activated dynamics to activated dynamics as 
the density increases and caging sets in. In this latter 
respect, the behavior of these systems is distinctly dif¬ 
ferent from that of a hard-sphere crystal with defects, 
where the particles are always caged and the diffusion of 
defects will occur via an activated process. 

There are, of course, some important differences. In 
the channel system, the growing bond-orientational or¬ 
der that we associate with the caging of disks is clearly 
visible in the structure factor, and this is not the case 
for the three-dimensional glass-forming liquids, where 
higher-order correlation functions are required to reveal 
bond-orientational order. The difference here is due to 
the channel walls, which form part of the cages: their di¬ 
rection is fixed and so constrains the possible directions 
of the “bonds” between nearest-neighbor disks. This can 
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be seen from the top and bottom diagrams Fig. 2, where 
only three nearest-neighbor bond directions are possible 
at (j) ^ 0max- 

Another difference appears in the behavior of the dy¬ 
namical susceptibility which, for the channel sys¬ 

tem, does not tend to zero for t ^ oo. We note, 
however, that this is the expected behavior of Xii^) in 
defect-mediated models of glassy dynamics in one dimen¬ 
sion [28]. 

We believe that our work strongly supports the com¬ 
mon idea [13-17] that glass behavior is a consequence 
of geometry and the local arrangements around the 
molecules in the supercooled liquid. Our system is 
sufficiently simple that we can quantitatively relate its 
dynamical features to structural features. The three- 
dimensional problem is much richer and success along 
these lines is probably only just starting [17]. 

The caging effects in our system mimic those seen in 
three dimensions. These are normally modeled by mode¬ 
coupling theory [31] and are associated with a genuine 
dynamical transition, within that approximation. Our 
system is effectively one dimensional and is unlikely to 
have a genuine phase transition. We are skeptical that 
the features that we see in the dynamics could be ex¬ 
plained in any way by mode-coupling arguments. They 


seem instead to be more naturally explained by dynam¬ 
ical processes associated with the developing structural 
order in the system. We suspect that the same might be 
true of three-dimensional systems. 

A noteworthy feature of the dynamics of our system of 
disks is the existence of two long time scales, r and 
In three dimensions, only the analogue of r, which is the 
cage-breaking time Tq,, is normally discussed. But in two 
dimensions a second, much longer, time scale has been 
revealed in experimental studies of polydisperse colloidal 
crystals, reported by Tanaka and co-workers [32, 33]. In 
their systems there is a growing lengthscale ^ for bond- 
orientational order, and their second long timescale rj is 
associated with dynamical correlations on that length- 
scale. This is strikingly similar to what we have found 
for the case of disks in a channel, where the growing 
lengthscale is that of zigzag ordering, which is a form of 
bond-orientational order. It might, of course, be that the 
second long time scale is a feature of the dynamics only 
for one- and two-dimensional systems. One would expect 
two long timescales to exist in three dimensions if escap¬ 
ing the cage could be associated with moving a defect; 
the longer time scale would then be associated with the 
time that the system needs to reach equilibrium via the 
motion of defects. 
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